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The Earth's core consists of a soUd bah with a radius of 1221 Km, surrounded by a Uquid shell 
which extends up to 3480 Km from the centre of the planet, roughly half way towards the surface 
(the mean radius of the Earth is 6373 km). The main constituent of the core is iron, and therefore 
the melting temperature of iron at the pressure encountered at the boundary between the solid 
and the liquid (the ICB) provides an estimate of the temperature of the core. Here I report the 
melting temperature of Fe at pressures near that of the ICB, obtained with first principles techniques 
based on density functional theory. The calculations have been performed by directly simulating 
solid and liquid iron in coexistence, and show that and at a pressure of ^ 328 GPa iron melts at 
^ 6370 ± 100 K. These findings are in good agreement with earlier simulations, which used exactly 
the same quantum mechanics techniques, but obtained melting properties from the calculation of 
the free energies of solid and liquid Fe [l] [2] [3j . 

PACS numbers: 



The study of iron under extreme conditions has a 
long hy story. In particular, numerous attempts have 
been made to obtain its high pressure melting proper- 
ties glElEllTllliaiTOllIIlIIl]. Experimentally, Earth's 
core conditions can only be reproduced by shock wave 
(SW) experiments, in which a high speed projectile is 
fired at an iron sample, and upon impact high pressure 
and high temperature conditions are produced. By vary- 
ing the speed of the projectile it is possible to investi- 
gate a characteristic pressure-volume relation known as 
the Hugoniot |13 , and even infer temperatures, although 
a word of caution here is in order, as temperature es- 
timates are often based on the knowledge of quantities 
like the constant volume specific heat and the Griineisen 
parameter, which are only approximately known at the 
relevant conditions [10]. If the speed of the projectile 
is high enough, the conditions of pressure and temper- 
ature are such that the sample melts, and it is there- 
fore possible to obtain points on the melting curve, of 
course with the caveat mentioned above about temper- 
ature measurements. An alternative route to high pres- 
sure high temperature properties is the use of diamond 
anvil cells (DAG), in which the sample is surrounded by a 
pressure medium and statically compressed between two 
diamond anvils. In DAG experiments pressure and tem- 
peratures can be directly measured, and therefore these 
techniques should in principle be more reliable to inves- 
tigate melting proeprties. Unfortunately, in the case of 
iron is it not so, and there is a fairly large range of results 
obtained by different groups HISl El H [3 19] . 

An alternative approach used for the past ten years or 
so has been to employ theory -and in particular quan- 
tum mechanics techniques based on density functional 
theory- to calculate the high pressure melting curve of 
iron. A number of groups have used different approaches 



to the problem. Our own strategy has been to calculate 
the Gibbs free energy of solid and liquid iron, and then 
obtain the melting curve by imposing their equality for 
any fixed pressure. We obtained a melting temperature 
of - 6350 K at 330 GPa [2]. The approach of Belonoshko 
et al. [14 was to fit an embedded atom model (EAM) to 
first principles calculations, and then calculate the melt- 
ing curve of the EAM. They obtained a temperature of 
~ 7050 K at 330 GPa. The approach of Laio et al. [15] 
was similar, altough they refitted their optimised model 
potential (0PM) to first principles calculation in a self- 
consistent way. They obtained a melting temperature of 
~ 5400 K at 330 GPa. We later re-conciled the results 
of Belonoshko et al [14^ with ours, by showing that the 
difference was due to a difference in free energies between 
their EAM and our DFT [16 . A similar argument would 
be responsible for the difference between our results and 
those of Laio et al [15^ . 

Here I am using an approach to melting which is 
independent from the free energy technique used ear- 
lier [H El [3], and the main motivation of this work is 
to provide an alternative route to the calculation of the 
melting properties of Fe. The method employed here is 
that of the coexistence of phases, in which solid and liq- 
uid iron are simulated in coexistence. The first time that 
the method was used in the context of first principles cal- 
culations was for the low pressure melting curve of alu- 
minium [17], where it was shown to deliver the same re- 
sults as the free energy method [18]. It was later applied 
to compute the melting curve of LiH [19] , hydrogen [20] 
and MgO [21 . 

The coexistence method is intrinsically expensive, as 
it requires large simulation cells and long simulations. It 
can be applied in a number of diffent ways, here I have 
used the NVE ensemble, i.e. constant number of atoms 
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A/", constant volume V and constant internal energy E. 
In the NVE ensemble, for each chosen volume V there 
is a whole range of energies E for which solid and liquid 
can coexist for long time; the average temperature and 
pressure along the simulation then provide a point on the 
melting curve. If the energy E is above(below) the range 
for which coexistence can be maintained, the system will 
completely melt (solidify), and the simulation does not 
provide useful melting properties informations. It should 
be pointed out that any finite system will eventually melt 
or solidify if simulated for long enough, due to sponta- 
neous fluctuations. However, melting(solidification) re- 
sulting from a too high (low) value of E typically appear 
on much shorter time scales. 

The present calculations have beed performed with 
density functional theory with the generalised gradient 
approximation known as PW91 [22] and the projector 
augmented wave method [24l [25] as implemented in the 
VASP code [23]. An efficient extrapolation of the charge 
density was employed ^7]. Single particle orbitals were 
expanded in plane waves with a cutoff of 300 eV, and I 
used the finite temperature implementation of DFT as 
developed by Mermin [28]. These settings are exactly 
equivalent to those used in our previous work [U (U [3^ , 
so the melting properties obtained here will be directly 
comparable to those early ones. The simulations have 
been performed on hexagonal closed packed (hep) cells 
containing 980 atoms (7 x 7 x 10), using the T point only. 
For the temperatures of interest here the use of the T 
point provides completely converged results. The time 
step in the molecular dynamics simulations was 1 femto- 
second (fs), and the self-consistency on the total energy 
2 X 10~^ eV. With these prescriptions the drift in the 
constant of motion was ~ 0.5 K/pico-second. 

The coexistence simulations were prepared by start- 
ing from a perfect hep crystal, which was initially ther- 
malised to ~ 6300 K for 1 pico-second (ps). Then half 
of the atoms in the cell were clamped and the tempera- 
ture was raised to a very high value, to melt the other 
half of the cell. Once a good melt was obtained, the 
temperature was reduced back to 6300 K and the system 
thermalised for one additional ps, after which the simu- 
lation was stopped, new initial velocities were assigned 
to the atoms and the simulation continued in the micro- 
canonical ensemble. The simulations were monitored us- 
ing the density profile, calculated by dividing the simula- 
tion cell in 100 slices parallel to the solid-liquid interface 
and counting the number of atoms in each slice; in the 
solid region this is a periodic function, with large num- 
ber of atoms if the slice coincide with an atomic plane, 
and small values if it falls between atomic layers. In the 
liquid region it fluctuates randomly around some average 
value. 

I performed five different simulations, starting with 
different amounts of internal energies E^ provided to 
the system by assigning different initial velocities to the 




FIG. 1: (Color online) Snapshot of a DFT molecular dynamics 
simulation showing solid and liquid iron in coexistence. The 
simulation cell contains 980 atoms. 
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FIG. 2: Temperature (upper panel) and pressure (lower 
panel) for a simulation of solid and liquid iron in coexistence. 

atoms. The simulation with the highest value of E com- 
pletely melted after ~ 6 ps. The one with the lowest 
amount of E solidified after ~ 11 ps. Among the other 
three, one melted after ~ 14 ps, one after ~ 24 ps, while 
the last one has remained in coexistence for the whole 
length of ~ 25 ps. However, most of these simulations 
were coexisting for long enough, so that useful melting 
information from the period of coexistence could actually 
be extracted in almost all cases. 

A snapshot of a simulation with solid and liquid in 
coexistence is show in Fig. [l] [29] . 

In Fig. |2] I display the temperatures and the pressures 
corresponding to the simulation that remained in coexis- 
tence for the whole 25 ps length, which provides a melting 
point (p,T) = (328 ±1 GPa, 6370 ±100 K). It is interest- 
ing to notice a temperature excursion in the simulation 
after ~ 15 ps, which lasts for ~ 5 ps. This tempera- 
ture variation is anti-correlated to a pressure variation, 
and corresponds to a temporary loss of some liquid in 
the cell, with latent heat of fusion converted into kinetic 
energy, and volume of fusion responsible for the drop in 
pressure. Large excursions of these type may provoke ac- 
cidental melting (or freezing), even if the internal energy 
E is within the range of coexistence. This problem is 
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FIG. 3: Temperature (upper panel) and pressure (lower 
panel) for a simulation of solid and liquid iron in coexistence. 
The system eventually completely solidifies, with a drop in 
pressure and an increase in temperature due to release of vol- 
ume and latent heat of fusion respectively. 
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FIG. 4: As in Fig. |3] but here the system eventually melts. 

mitigated by the use of large simulation cells, and there- 
fore this is one of the reasons why large simulation cells 
are needed in conjunction with the coexistence approach. 

In Fig. [3] I show a simulation that eventually solidi- 
fied, however, as mentioned above, coexistence was main- 
tained for a long period, and the information gathered 
by the central part of the simulation can still be used 
to obtain a point on the melting curve, and the result is 
(p, T) = (324 ± 1 GPa, 6250 ± 100 K), which is consistent 
with the previous point. 

Similarly, in Fig. [4] I show one of the simulations that 
eventually melted, and by taking the average tempera- 
ture and pressure from the central part of the simulation 
I get (p, T) = (331 ± 1 GPa, 6430 ± 100 K), which is also 
consistent with the other previous two points. 

The c/a value used in the simulations was fixed at 
1.6, which resulted in slightly non- hydrostatic condi- 
tions. To study the effect of non-hydrostaticity, and 
that of the relatively small size of the simulation cells. 



I have performed simulations using a classical embed- 
ded atom model (EAM), adapted to deliver results very 
close to the present ab-inition techniques [16 . I have 
performed simulations on cells containing 7840 atoms 
(14 X 14 X 20), and found that at a pressure of 324 GPa 
the effect of using a small cell containing only 980 atoms 
is to raise the melting temperature by ~ 100 K. The 
ab-initio non-hydrostatic conditions with c/a = 1.6 are 
similarly reproduced by the EAM, which also shows that 
with c/a = 1.65 the simulations are almost exactly un- 
der hydrostatic conditions [30 . The effect of the non- 
hydrostaticity with c/a = 1.6 is to reduce the melting 
temperature by ~ 100 K, so that the combined effects of 
non-hydrostaticity and small size cancel each other. A 
final check was performed by repeating the simulations 
using 62720 atom-cells (28 x 28 x 40), which showed es- 
sentially no differences with the results obtained using 
7840 atom-cells. 

All the present first-principles coexistence results are 
displayed in Fig. [5) which also contains experimental 
and previous theoretical results. The filled square cor- 
responds to the simulation which maintained coexistence 
throughout (Fig.|2|, while the empty squares to the other 
4 simulations, for which the final part has been discarded. 
As noted above, it is clear that all simulations provide 
similar melting points, with the exception perhaps of the 
point correponding to the simulation that melted only af- 
ter ~ 6 ps (point at highest temperature in Fig[5|. The 
comparison of the present results with the earlier melting 
curve obtained using the free energy approach [2] shows 
excellent agreement, and the two sets of data therefore 
support each other. I also report on the same figure the 
DAG experiments of Refs. 01 HI H [71 [9], the SW ex- 
periments of Refs. [lOl [m [12] and the calculations of 
Refs. |T5]. As mentioned above, we identified the 
reasons of the differences between our melting curve and 
that calculated by Belonoshko et al. [14] being the free en- 
ergy differences between their EAM and our DFT. Once 
these differences are taken into account it is possible to 
"correct" the EAM melting curve. The two red dots on 
the figure show the corrected EAM results at two differ- 
ent pressures, which agree with our DFT melting curve. 

In conclusion, I have presented here calculations for 
the melting temperature of iron under Earth's core con- 
ditions, obtained with density functional theory and the 
technique of the coexistence of phases. The DFT tech- 
niques are the same as those employed in our earlier cal- 
culations, where we computed the melting curve using 
the free energy approach [2] |3]. The consistence be- 
tween the present calculations and those early ones is 
expected, because the DFT technicalities are the same, 
and provides an independent check on the accuracy of 
those early free energy calculations. The DFT melting 
temperature of Fe at a pressure of ^ 330 GPa is therefore 
confirmed to be in the region of ~ 6300 — 6400 K. 
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FIG. 5: (Color online) Comparison of melting curve of Fe 
from present calculations with previous experimental and ab 
initio results; blue filled square: DFT coexisting simulation 
for over 25 ps; blue open squares: DFT coexisting simulations 
ending up to complete solids or liquids, results are gathered 
from the region of coexistence; black heavy solid curve: DFT 
melting curve from Ref. |2j long-dashed blue curve: DFT- 
0PM results of Ref. [E]; hght purple sohd curve: DFT-EAM 
results of Ref. [14]; filled red circles: DFT-EAM "corrected" 
results (see text); black chained and maroon dashed curves: 
DAC measurements of Refs. [4] and [5]; green open diamonds: 
DAC measurements of Ref. 6^; green plus: DAC measurement 
of Ref. 9 ; green filled triangle: DAC measurement of Ref. [7 ; 
black stars, black open circle and pink filled diamond: shock 
experiments of Refs. [12], and \1X». Error bars are those 
quoted in original references. 



The next question that one might ask is how accu- 
rate is DFT for this problem. We argued in our pre- 
vious work that DFT should indeed be quite reliable, 
mainly because solid and liquid iron have very similar 
structures, so that possible errors would largely cancel 
between the two phases. However, we also pointed out 
that DFT does not seem to reproduce the zero tempera- 
ture pressure- volume equation of state of hep iron com- 
pletely correctly, possibly underestimating the pressure 
by ~ 2.5%. We then argued that this error could propa- 
gate to the melting curve, resulting in a lowering of tem- 
peratures which at a pressure of 330 GPa could be in 
the region of ~ 150 K [2 . This would bring the melting 
temperature of Fe at ICB condition to ~ 6200 K. It will 
be interesting to re-visit this problem with more accurate 
quantum mechanics techniques, and we are planning to 
do so by using quantum Monte Carlo. We will report on 
these results in due course. 

This work was conducted as part of a EURYI scheme 
award as provided by EPSRC (see www.esf.org/euryi). 
Calculations have been performed on the UK national 
facility HECToR, using allocation of time from the Min- 
eral Consortium and from a EPSRC Capability Chal- 



lenge grant. Calculations were also performed on the 
UCL research-computing facility Legion, and initially on 
the Cambridge High Performance facility Darwin. Sim- 
ulations were typically run on 256 cores, each molecular 
dynamics step of 1 fs taking ~ 7.5 minutes. 
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